Independent Research Project: Fitting Sersic Profiles to Elliptical and Spiral Galaxy Observations using SEO

Jack Colvin, Along with my group, Larrence Xing, Kevin Wu, Annabelle Yarborough and Olivia Lucal

1. Introduction

Sersic Profiles, Introduced by Argentinian astronomer José Luis Sérsic, are a generalization of de Vacouleurs' law, which describes the light prodiles of elliptical galaxies, to all galaxy types. It specically fits galaxies to a symmetrical light distribution described by the equation $I * e^{-b((\frac{r}{r_{e}})^{\frac{1}{n}}-1)}$, where n is a varible index that describes how quickly the light in a given galaxy distribution falls off from its max value, with a higher index representing a sharper fall-off. This is useful in studying galaxies specifically in how spiral and elliptical galaxies are sharply dilineated with respect to this index: eliptical galaxies generally have sersic index of around 4, while spiral galaxies have an index of around 1 (which is an exponential profile). With this in mind, I took 6 observations, 3 of spiral galaxies (m51, m81, m101) and three eliptical gal;axies (m89, m49, m60) to see if I could

1. Retrieve an index of ~4 for elliptical galaxies and ~1 for spiral galaxies
2. And if I could not, whether I could at least measure a difference between the two galaxy types, where eliptical galaxies have higher sersic indices and thus a more concentrated light distribution, and spirals have a lower index and thus a more diffuse light distribution

2. Methods and Data Collection

I started by taking 120s exposures of each galaxy with seo, which represented a good balance betweenb having a high enough SNR, while avoiding the pixels within the galaxies from spilling into neigboring pixels (which would cause big problems whgen trying to conduct sersic fits). I took exposures in all three bands in order to be able to create pretty rgb images for all three galaxies! Finaly, I constructed a sersic function and a least sqaures residual function in order to use least sqaures fitting to fit a sersic model to each galaxy.

While this relatively simple model worked for elliptical galaxies relatively well, I had to make some modifications when working with spiral galaxies. First and foremost, the bulge is both much brighter than the rest of the disk, and has a different, more concentrated sersic profile, which was resulting in a one-component model only fitting to the bulge, and not to the rest of the galaxy, and thus overestimating the overall sersic index. To correct for this, I switched to a two-component model when dealing with disks. I also suspect that the spiral arms are also caussing fitting issues, as they are assymetric, and are likely thus ignored by the model. I hav3 enot yet corrected for this effect, but I will certainly do so as I continue to work on this project in the future. Finally, I tried masking out the bulge entirely, to see if doiung so would give me a better or worse result than the two component model, with varying results, and I also tried fitting a two component model with a sersic index for the bulge set between 3-6 as another method to remove the variability of the bulge as my goal is to fit a model to the disk halo.

import numpy as np
from astropy.io import fits
import matplotlib.pyplot as plt
from astropy.visualization import ZScaleInterval
from matplotlib.colors import LogNorm
from astropy.nddata import Cutout2D
from astropy.modeling import models, fitting 
from astropy.visualization.lupton_rgb import make_lupton_rgb
from astropy.visualization import LogStretch
from scipy.optimize import least_squares
import matplotlib.pyplot as plt

def plot_prettier(dpi=150, fontsize=11, usetex=False): 
    '''
    Make plots look nicer compared to Matplotlib defaults
    Parameters: 
        dpi - int, "dots per inch" - controls resolution of PNG images that are produced
                by Matplotlib
        fontsize - int, font size to use overall
        usetex - bool, whether to use LaTeX to render fonds of axes labels 
                use False if you don't have LaTeX installed on your system
    '''
    plt.rcParams['figure.dpi']= dpi
    plt.rc("savefig", dpi=dpi)
    plt.rc('font', size=fontsize)
    plt.rc('xtick', direction='in') 
    plt.rc('ytick', direction='in')
    plt.rc('xtick.major', pad=5) 
    plt.rc('xtick.minor', pad=5)
    plt.rc('ytick.major', pad=5) 
    plt.rc('ytick.minor', pad=5)
    plt.rc('lines', dotted_pattern = [2., 2.])
    if usetex:
        plt.rc('text', usetex=usetex)
    else:
        plt.rcParams['mathtext.fontset'] = 'cm'
        plt.rcParams['font.family'] = 'serif'
        plt.rcParams['font.serif'] = ['Times New Roman'] + plt.rcParams['font.serif']
        

plot_prettier(dpi=150, fontsize=11)

Model Fitting for M51

iband = fits.open("/Users/jackcolvin//m51_i-band_120.0s_bin1_250603_054315_jackcolvin_seo_0_FCAL.fits")
gband = fits.open("/Users/jackcolvin//m51_g-band_120.0s_bin1_250603_054042_jackcolvin_seo_0_FCAL.fits")
rband = fits.open("/Users/jackcolvin//m51_r-band_120.0s_bin1_250603_053809_jackcolvin_seo_0_FCAL.fits")

r = iband[0].data
g = rband[0].data
b = gband[0].data

#lupton_rgb doesn't support scaling
#so that must be done manually
zscale = ZScaleInterval()
r_scaled = zscale(r)
g_scaled = zscale(g)
b_scaled = zscale(b)



image = make_lupton_rgb(r_scaled, g_scaled, b_scaled, Q = 1, stretch = .9, minimum = .3)
plt.imshow(image, origin = 'lower')
plt.show()

#Now we must find a best fit sersic profile
#and see if our best fit index fits our prediction 
#based on whether this galaxy is spiral or elliptical. 
#also check to see how different the function is across different bands
from scipy.optimize import least_squares
import matplotlib.pyplot as plt

scale_factor = 10000 #rescale calibrated fluxes as they are so small
cutout = Cutout2D(r, (1100, 1200), (1000, 1000))
data0 = cutout.data
data = (data0-np.median(r))*scale_factor
interval = ZScaleInterval()
vmin,vmax = interval.get_limits(data)
plt.imshow(data, origin = 'lower', cmap = 'gray', vmin = vmin, vmax = vmax )


y = np.zeros(shape = data.shape, dtype = int)
x = np.zeros(shape = data.shape, dtype = int)
for i in range(len(y[0])):
    y[i] += len(y[0]) - (i+1)
    x[:,i] += i

    
#define an example sersic profile:
def sersic_2d(x, y, I_e, R_e, n, x0, y0):
    r = np.sqrt((x - x0)**2 + (y - y0)**2)
    b_n = 2 * n - 1/3 + (4/(405*n)) + (46/(25515*n**2))
    return I_e * np.exp(-b_n * ((r / R_e)**(1/n) - 1))

def residuals(params, x, y, data):
    I_e, R_e, n, x0, y0 = params
    model = sersic_2d(x, y, I_e, R_e, n, x0, y0)
    return (model - data).ravel()

p0 = [np.max(data), 50, 1.0, 500, 500] #guess the expected sersic index for a spiral galaxy
result = least_squares(residuals, p0, args=(x, y, data), bounds=([0, 1, 0.1, 0, 0], [np.inf, np.inf, 10, 500, 500]))

I_e_fit, R_e_fit, n_fit, x0_fit, y0_fit = result.x
print(f"Fit: I_e = {I_e_fit:.5f}, R_e = {R_e_fit:.2f}, n = {n_fit:.2f}, x0 = {x0_fit:.2f}, y0 = {y0_fit:.2f}")

Fit: I_e = 0.04911, R_e = 194.49, n = 2.33, x0 = 493.56, y0 = 500.00

model_image = sersic_2d(x, y, *result.x)

import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno',  vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm',  vmin = vmin, vmax = vmax )
plt.show()

#we can now evaluate the chi2 of our fit: 
data1 = data.flatten()
mask = data1 > 0
data2 = data1[mask]
model_image1 = model_image.flatten()
model_image2 = model_image1[mask]
residual = (data2 - model_image2)
uncertainty = (np.sqrt(data2))#use poisson uncertainty on fluxes
chi2 = np.sum((residual/uncertainty)**2)
dof = len(data2) - len(result.x)
chi2_red = chi2/ dof
print(f"Chi-squared: {chi2:.2f}")
print(f"Reduced Chi-squared: {chi2_red:.2f}")

#very low reduced chi2, almost symptomatic of over fitting 
#also, our sersic index is very high, almost 4
#which indicates that it is capturing only the bulge, rather than the full disk
#so now I wanted to try a two-component model

Chi-squared: 59764.06
Reduced Chi-squared: 0.10

def two_component_model(x, y, params):
    # params = [I_e_bulge, R_e_bulge, n_bulge, I_e_disk, R_e_disk, x0, y0]
    I_e_bulge, R_e_bulge, n_bulge, I_e_disk, R_e_disk, n_disk, x0, y0 = params
    bulge = sersic_2d(x, y, I_e_bulge, R_e_bulge, n_bulge, x0, y0)
    disk = sersic_2d(x, y, I_e_disk, R_e_disk, n_disk, x0, y0)  
    return bulge + disk

def residuals_two_components(params, x, y, data):
    model = two_component_model(x, y, params)
    return (model - data).ravel()

p1 = [np.max(data), 50, 4, np.max(data)/1.5, 200, 1, 500, 500]
bounds_lower = [0, 1, .5, 0, 10, .5, 0, 0]
bounds_upper = [np.inf, 200, 8, np.inf, 500, 8, data.shape[1], data.shape[0]]


result = least_squares(residuals_two_components, p1, args=(x, y, data), bounds=(bounds_lower, bounds_upper))

I_e_bulge_fit, R_e_bulge_fit, n_bulge_fit, I_e_disk_fit, R_e_disk_fit, n_disc_fit, x0_fit, y0_fit = result.x

print(f"Fit results:")
print(f"Bulge: I_e={I_e_bulge_fit:.2f}, R_e={R_e_bulge_fit:.2f}, n={n_bulge_fit:.2f}")
print(f"Disk: I_e={I_e_disk_fit:.2f}, R_e={R_e_disk_fit:.2f}, n={n_disc_fit:.2f}")
print(f"Center: x0={x0_fit:.2f}, y0={y0_fit:.2f}")

Fit results:
Bulge: I_e=0.72, R_e=18.74, n=1.05
Disk: I_e=0.07, R_e=185.16, n=0.77
Center: x0=494.71, y0=515.74

model_image = two_component_model(x, y, np.array(result.x))

import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm', vmin = vmin, vmax = vmax)
plt.show()

#xlearly, the sersic profile is beinmg affected by the spiral arms
#which is shown in the residuals plot
#thus masking pout the spiral arms may be a necessary step to getting a better fit. 

#now try fitting to the disk profile by removing the center bulge entirely 
#to see if the fit becomes becomes better after doing so

from astropy.convolution import Gaussian2DKernel, convolve

kernel = Gaussian2DKernel(x_stddev=20)  # tune stddev to smooth scale
smooth = convolve(data, kernel)


residual = data - smooth

sigma = np.std(residual)
threshold = sigma

R_bulge = 50  # adjust based on your galaxy
yy, xx = np.indices(data.shape)
r = np.sqrt((xx - 500)**2 + (yy - 500)**2)
central_mask = r < R_bulge

# Mask = True where pixels are *not* spiral arms (we want to keep them)
#but we also need to avoid masking out the central bulge. 
mask0 = np.abs(residual) < threshold
mask = mask0


data_corr = np.where(mask, data, np.median(data))


y = np.zeros(shape = data_corr.shape, dtype = int)
x = np.zeros(shape = data_corr.shape, dtype = int)
for i in range(len(y[0])):
    y[i] += len(y[0]) - (i+1)
    x[:,i] += i



p1 = [np.max(data), 50, 4, np.max(data)/1.5, 200, 1, 500, 500]
bounds_lower = [0, 1, .5, 0, 10, .5, 0, 0]
bounds_upper = [np.inf, 200, 8, np.inf, 500, 8, data.shape[1], data.shape[0]]


result = least_squares(residuals_two_components, p1, args=(x, y, data_corr), bounds=(bounds_lower, bounds_upper))

I_e_bulge_fit, R_e_bulge_fit, n_bulge_fit, I_e_disk_fit, R_e_disk_fit, n_disc_fit, x0_fit, y0_fit = result.x

print(f"Smoothed Fit results:")
print(f"Bulge: I_e={I_e_bulge_fit:.2f}, R_e={R_e_bulge_fit:.2f}, n={n_bulge_fit:.2f}")
print(f"Disk: I_e={I_e_disk_fit:.2f}, R_e={R_e_disk_fit:.2f}, n={n_disc_fit:.2f}")
print(f"Center: x0={x0_fit:.2f}, y0={y0_fit:.2f}")

Smoothed Fit results:
Bulge: I_e=0.00, R_e=1.00, n=0.50
Disk: I_e=0.06, R_e=178.04, n=0.84
Center: x0=486.89, y0=529.90

model_image = two_component_model(x, y, np.array(result.x))

import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Smooth_data')
plt.imshow(data_corr, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data_corr - model_image, origin='lower', cmap='coolwarm', vmin = vmin, vmax = vmax)
plt.show()

#xlearly, the sersic profile is beinmg affected by the spiral arms
#which is shown in the residuals plot
#While the bulge plot is no longer accurate 
#as the bulge is filtered out 
#the model sersic fit is much closer to what we expect 
#and the residual plot looks much bteter for the outside. 

iband.close()
rband.close()
gband.close()

Basic Analysis for M51

• Our first model fits just for the bulge, and retrieves a sersic index of 2.33, which is fairly accurate for the bulge of a spiral galaxy, but obviously a severe over-estimation for the disc halo. This can be seen in the model plot, which fits very well for the bulge, but doesn't consider the light halo of the rest of the galaxy at all.
• We get a slightly better estimation when switching to a 2 component model, with a value of .77 for the disc halo, which is an underestimation, but much more reasonable. The two component model, however, also severly underestimates the index of the bulge, giving a value of only 1.05, which is far too diffuse.
• The best value may actually be from a combinatiomn of these two models, which could be an interesting hypothesis to test in an extension to this project.
• Finally, the best value we retrieve for the disc is when we remove the vbulge entirely, along with smoothing away the spiral arms slightly, which gives us a value of .84, which is very reasonable, especially considering that SEO has very low angular resoltion,a nd thus it is reasonable to assume that some smoothing and spreading out of the light profile may be taking place. This fit also has the best residual plot, apart form the center, suggesting it may be the most accurate of the three methods, at least for fitting the disk model.
• All three models fail to capture the spiral arms in any meaningful way, as evidenced by the residual plots.

Model Fitting for M101

iband0 = fits.open("/Users/jackcolvin//m101_i-band_120.0s_bin1_250603_053348_jackcolvin_seo_0_FCAL.fits")
gband0 = fits.open("/Users/jackcolvin/m101_g-band_120.0s_bin1_250603_053050_jackcolvin_seo_0_FCAL.fits")
rband0 = fits.open("/Users/jackcolvin//m101_r-band_120.0s_bin1_250603_052816_jackcolvin_seo_0_FCAL.fits")

r0 = iband0[0].data
g0 = rband0[0].data
b0 = gband0[0].data

#lupton_rgb doesn't support scaling
#so that must be done manually
zscale = ZScaleInterval()
r_scaled0 = zscale(r0)
g_scaled0 = zscale(g0)
b_scaled0 = zscale(b0)



image = make_lupton_rgb(r_scaled0, .9*g_scaled0, b_scaled0, Q = 1, stretch = .7, minimum = .25)
plt.imshow(image, origin = 'lower')
plt.show()

#Now we must find a best fit sersic profile
#and see if our best fit index fits our prediction 
#based on whether this galaxy is spiral or elliptical. 
#also check to see how different the function is across different bands
from scipy.optimize import least_squares
import matplotlib.pyplot as plt

scale_factor = 10000 #rescale calibrated fluxes as they are so small
cutout = Cutout2D(r0, (1100, 1250), (600, 600))
data0 = cutout.data
data = (data0-np.median(r0))*scale_factor
interval = ZScaleInterval()
vmin,vmax = interval.get_limits(data)
plt.imshow(data, origin = 'lower', cmap = 'gray', vmin = vmin, vmax = vmax )


y = np.zeros(shape = data.shape, dtype = int)
x = np.zeros(shape = data.shape, dtype = int)
for i in range(len(y[0])):
    y[i] += len(y[0]) - (i+1)
    x[:,i] += i

    
#define an example sersic profile:
def sersic_2d(x, y, I_e, R_e, n, x0, y0):
    r = np.sqrt((x - x0)**2 + (y - y0)**2)
    b_n = 2 * n - 1/3 + (4/(405*n)) + (46/(25515*n**2))
    return I_e * np.exp(-b_n * ((r / R_e)**(1/n) - 1))

def residuals(params, x, y, data):
    I_e, R_e, n, x0, y0 = params
    model = sersic_2d(x, y, I_e, R_e, n, x0, y0)
    return (model - data).ravel()

p0 = [np.max(data), 50, 1.0, 250, 300] #guess the expected sersic index for a spiral galaxy
result = least_squares(residuals, p0, args=(x, y, data), bounds=([0, 1, 0.1, 0, 0], [np.inf, np.inf, 10, 500, 500]))

I_e_fit, R_e_fit, n_fit, x0_fit, y0_fit = result.x
print(f"Fit: I_e = {I_e_fit:.5f}, R_e = {R_e_fit:.2f}, n = {n_fit:.2f}, x0 = {x0_fit:.2f}, y0 = {y0_fit:.2f}")

Fit: I_e = 0.00152, R_e = 2192.04, n = 4.42, x0 = 237.36, y0 = 298.53

model_image = sersic_2d(x, y, *result.x)

import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno',  vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm',  vmin = vmin, vmax = vmax )
plt.show()

p1 = [np.max(data), 10, 4, np.max(data)/1.5, 200, 1, 300, 280]
bounds_lower = [0, 1, 2, 0, 10, .5, 0, 0]
bounds_upper = [np.inf, 200, 8, np.inf, 500, 8, data.shape[1], data.shape[0]]


result = least_squares(residuals_two_components, p1, args=(x, y, data), bounds=(bounds_lower, bounds_upper))

I_e_bulge_fit, R_e_bulge_fit, n_bulge_fit, I_e_disk_fit, R_e_disk_fit, n_disc_fit, x0_fit, y0_fit = result.x

print(f"Fit results:")
print(f"Bulge: I_e={I_e_bulge_fit:.2f}, R_e={R_e_bulge_fit:.2f}, n={n_bulge_fit:.2f}")
print(f"Disk: I_e={I_e_disk_fit:.2f}, R_e={R_e_disk_fit:.2f}, n={n_disc_fit:.2f}")
print(f"Center: x0={x0_fit:.2f}, y0={y0_fit:.2f}")

Fit results:
Bulge: I_e=0.06, R_e=48.55, n=2.00
Disk: I_e=0.04, R_e=270.66, n=0.92
Center: x0=237.33, y0=298.52

model_image = two_component_model(x, y, np.array(result.x))

import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm', vmin = vmin, vmax = vmax)
plt.show()

#xlearly, the sersic profile is beinmg affected by the spiral arms
#which is shown in the residuals plot
#thus masking pout the spiral arms may be a necessary step to getting a better fit. 

def two_component_ellipse(x, y, params):
    # params = [I_e_bulge, R_e_bulge, n_bulge, q_bulge, theta_bulge, I_e_disk, R_e_disk, q_disk, theta_disk, x0, y0]
    I_e_bulge, R_e_bulge, n_bulge, q_bulge, theta_bulge, I_e_disk, R_e_disk, n_disk, q_disk, theta_disk, x0, y0 = params
    bulge = sersic_2d_elliptical(x, y, I_e_bulge, R_e_bulge, n_bulge, x0, y0, q_bulge, theta_bulge)
    disk = sersic_2d_elliptical(x, y, I_e_disk, R_e_disk, n_disk, x0, y0, q_disk, theta_disk)  
    return bulge + disk

def residuals_two_ellipse(params, x, y, data):
    model = two_component_ellipse(x, y, params)
    return (model - data).ravel()


p1 = [np.max(data), 15, 4, .5, 1, np.max(data)/2, 200, 1, .5, 1, 230, 270]
bounds_lower = [0, 0, 2, 0, 0, 0, 50, .5, 0, 0, 0, 0]
bounds_upper = [np.inf, 50, 6, 1, 2*np.pi, np.inf, 400, 8, 1, 2*np.pi, data.shape[1], data.shape[0]]


result = least_squares(residuals_two_ellipse, p1, args=(x, y, data), bounds=(bounds_lower, bounds_upper))

I_e_bulge_fit, R_e_bulge_fit, n_bulge_fit, q_bulge_fit, theta_bulge_fit, I_e_disk_fit, R_e_disk_fit, n_disc_fit, q_disc_fit, theta_disc_fit, x0_fit, y0_fit = result.x

print(f"Fit results:")
print(f"Bulge: I_e={I_e_bulge_fit:.2f}, R_e={R_e_bulge_fit:.2f}, n={n_bulge_fit:.2f}, q ={q_bulge_fit:.2f}, theta = {theta_bulge_fit:.2f}")
print(f"Disk: I_e={I_e_disk_fit:.2f}, R_e={R_e_disk_fit:.2f}, n={n_disc_fit:.2f}, q={q_disc_fit:.2f}, theta={theta_disc_fit:.2f}")
print(f"Center: x0={x0_fit:.2f}, y0={y0_fit:.2f}")

Fit results:
Bulge: I_e=0.07, R_e=50.00, n=2.00, q =0.81, theta = 5.98
Disk: I_e=0.03, R_e=313.05, n=1.01, q=0.80, theta=2.54
Center: x0=237.30, y0=298.50

model_image = two_component_ellipse(x, y, np.array(result.x))

import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm', vmin = vmin, vmax = vmax)
plt.show()

#xlearly, the sersic profile is beinmg affected by the spiral arms
#which is shown in the residuals plot
#thus masking pout the spiral arms may be a necessary step to getting a better fit. 

iband0.close()
gband0.close()
rband0.close()

Basic Analysis for M101

• Our first model fits just fpor the bulge, and retrieves a sersic index of 4.42, which again is a reasonable estimate for the bulge, especially considering that m101 has a more compact central bulge, but again represents a poor model for the overall galaxy
• We get a very accurate model when seperating the bulge and the dic into a two component model, with an estimated index of .92 for the disc halo, which is very close to the expected value of .92.
• Yet again, however, the index for the bulge is underestimated, with a calculated value of n=2, which, while closer to being accurate than the value obtained for m51, is still quite far away from the expected value of 3-6.
• The last fit i tried was fixing the bulge to be > 2, which still gave me an underestimated index for my bulge, but gave me an index of 1.01 for my disc, which is extremely accurate.
• Finally, yet again, the residuals plot indicates that the fit is quite accurate, but does not take into account the spiral arms due to their lack of symmetry.

Model Fitting for M81

iband01 = fits.open("/Users/jackcolvin//m81_i-band_120.0s_bin1_250603_055612_jackcolvin_seo_0_FCAL.fits")
gband01 = fits.open("/Users/jackcolvin//m81_g-band_120.0s_bin1_250603_055340_jackcolvin_seo_0_FCAL.fits")
rband01 = fits.open("/Users/jackcolvin//m81_r-band_120.0s_bin1_250603_055109_jackcolvin_seo_0_FCAL.fits")

r01 = iband01[0].data
g01 = rband01[0].data
b01 = gband01[0].data

#lupton_rgb doesn't support scaling
#so that must be done manually
zscale = ZScaleInterval()
r_scaled01 = zscale(r01)
g_scaled01 = zscale(g01)
b_scaled01 = zscale(b01)



image = make_lupton_rgb(r_scaled01, .85*g_scaled01, b_scaled01, Q = 1, stretch = .7, minimum = .25)
plt.imshow(image, origin = 'lower')
plt.show()

scale_factor = 10000 #rescale calibrated fluxes as they are so small
cutout = Cutout2D(r01, (1000, 1200), (1000, 1000))
data0 = cutout.data
data = (data0-np.median(r01))*scale_factor
interval = ZScaleInterval()
vmin,vmax = interval.get_limits(data)
plt.imshow(data, origin = 'lower', cmap = 'gray', vmin = vmin, vmax = vmax )


y = np.zeros(shape = data.shape, dtype = int)
x = np.zeros(shape = data.shape, dtype = int)
for i in range(len(y[0])):
    y[i] += len(y[0]) - (i+1)
    x[:,i] += i

    
#define an example sersic profile:
def sersic_2d(x, y, I_e, R_e, n, x0, y0):
    r = np.sqrt((x - x0)**2 + (y - y0)**2)
    b_n = 2 * n - 1/3 + (4/(405*n)) + (46/(25515*n**2))
    return I_e * np.exp(-b_n * ((r / R_e)**(1/n) - 1))

def residuals(params, x, y, data):
    I_e, R_e, n, x0, y0 = params
    model = sersic_2d(x, y, I_e, R_e, n, x0, y0)
    return (model - data).ravel()

p0 = [np.max(data), 50, 1.0, 500, 500] #guess the expected sersic index for a spiral galaxy
result = least_squares(residuals, p0, args=(x, y, data), bounds=([0, 1, 0.1, 0, 0], [np.inf, np.inf, 10, 600, 600]))

I_e_fit, R_e_fit, n_fit, x0_fit, y0_fit = result.x
print(f"Fit: I_e = {I_e_fit:.5f}, R_e = {R_e_fit:.2f}, n = {n_fit:.2f}, x0 = {x0_fit:.2f}, y0 = {y0_fit:.2f}")

Fit: I_e = 0.20730, R_e = 148.97, n = 3.14, x0 = 501.55, y0 = 491.44

model_image = sersic_2d(x, y, *result.x)

import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno',  vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm',  vmin = vmin, vmax = vmax )
plt.show()

#a circularly symetric model does a bad job for this galaxy
#so now we must define an ellipse

def sersic_2d_elliptical(x, y, I_e, R_e, n, x0, y0, q, theta):
    x_shift = x - x0
    y_shift = y - y0
    x_rot = x_shift * np.cos(theta) + y_shift * np.sin(theta)
    y_rot = -x_shift * np.sin(theta) + y_shift * np.cos(theta)
    R = np.sqrt(x_rot**2 + (y_rot / q)**2)
    b_n = 2 * n - 1/3 + (4/(405*n)) + (46/(25515*n**2))
    return I_e * np.exp(-b_n * ((R / R_e)**(1/n) - 1))

def residuals_elliptical(params, x, y, data):
    I_e, R_e, n, x0, y0, q, theta = params
    model = sersic_2d_elliptical(x, y, I_e, R_e, n, x0, y0, q, theta)
    return (model - data).ravel()

p0 = [np.max(data), 50, 1.0, 500, 500, .5, np.pi/3] #guess the expected sersic index for a spiral galaxy
result = least_squares(residuals_elliptical, p0, args=(x, y, data), bounds=([0, 1, 0.1, 0, 0, 0, 0], [np.inf, np.inf, 10, 600, 600, 1, 2*np.pi]))

I_e_fit, R_e_fit, n_fit, x0_fit, y0_fit, q_fit, theta_fit = result.x
print(f"Fit: I_e = {I_e_fit:.5f}, R_e = {R_e_fit:.2f}, n = {n_fit:.2f}, x0 = {x0_fit:.2f}, y0 = {y0_fit:.2f}, q = {q_fit:.2f}, theta = {theta_fit:.2f}")

Fit: I_e = 0.20369, R_e = 181.48, n = 3.14, x0 = 501.58, y0 = 491.36, q = 0.71, theta = 0.99

model_image = sersic_2d_elliptical(x, y, *result.x)

import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno',  vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm',  vmin = vmin, vmax = vmax )
plt.show()

p1 = [np.max(data), 50, 4, .5, 1, np.max(data)/1.5, 200, 1, .5, 1, 500, 500]
bounds_lower = [0, 1, .5, 0, 0, 0, 10, .5, 0, 0, 0, 0]
bounds_upper = [np.inf, 200, 8, 1, 2*np.pi, np.inf, 500, 8, 1, 2*np.pi, data.shape[1], data.shape[0]]


result = least_squares(residuals_two_ellipse, p1, args=(x, y, data), bounds=(bounds_lower, bounds_upper))

I_e_bulge_fit, R_e_bulge_fit, n_bulge_fit, q_bulge_fit, theta_bulge_fit, I_e_disk_fit, R_e_disk_fit, n_disc_fit, q_disc_fit, theta_disc_fit, x0_fit, y0_fit = result.x

print(f"Fit results:")
print(f"Bulge: I_e={I_e_bulge_fit:.2f}, R_e={R_e_bulge_fit:.2f}, n={n_bulge_fit:.2f}, q ={q_bulge_fit:.2f}, theta = {theta_bulge_fit:.2f}")
print(f"Disk: I_e={I_e_disk_fit:.2f}, R_e={R_e_disk_fit:.2f}, n={n_disc_fit:.2f}, q={q_disc_fit:.2f}, theta={theta_disc_fit:.2f}")
print(f"Center: x0={x0_fit:.2f}, y0={y0_fit:.2f}")

Fit results:
Bulge: I_e=3.00, R_e=14.68, n=0.98, q =0.83, theta = 0.71
Disk: I_e=0.22, R_e=206.61, n=1.75, q=0.60, theta=1.06
Center: x0=501.58, y0=491.31

model_image = two_component_ellipse(x, y, np.array(result.x))

import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm', vmin = vmin, vmax = vmax)
plt.show()

p1 = [np.max(data), 50, 5, .5, 1, np.max(data)/2, 200, 1, .5, 1, 500, 500]
bounds_lower = [0, 10, 2.5, 0, 0, 0, 0, .5, 0, 0, 0, 0]
bounds_upper = [np.inf, 50, 6, 1, 2*np.pi, np.inf, 500, 8, 1, 2*np.pi, data.shape[1], data.shape[0]]


result = least_squares(residuals_two_ellipse, p1, args=(x, y, data), bounds=(bounds_lower, bounds_upper))

I_e_bulge_fit, R_e_bulge_fit, n_bulge_fit, q_bulge_fit, theta_bulge_fit, I_e_disk_fit, R_e_disk_fit, n_disc_fit, q_disc_fit, theta_disc_fit, x0_fit, y0_fit = result.x

print(f"Fit results:")
print(f"Bulge: I_e={I_e_bulge_fit:.2f}, R_e={R_e_bulge_fit:.2f}, n={n_bulge_fit:.2f}, q ={q_bulge_fit:.2f}, theta = {theta_bulge_fit:.2f}")
print(f"Disk: I_e={I_e_disk_fit:.2f}, R_e={R_e_disk_fit:.2f}, n={n_disc_fit:.2f}, q={q_disc_fit:.2f}, theta={theta_disc_fit:.2f}")
print(f"Center: x0={x0_fit:.2f}, y0={y0_fit:.2f}")

Fit results:
Bulge: I_e=1.33, R_e=10.00, n=2.50, q =0.02, theta = 3.22
Disk: I_e=0.20, R_e=182.42, n=3.14, q=0.71, theta=1.00
Center: x0=501.57, y0=491.45

model_image = two_component_ellipse(x, y, np.array(result.x))

import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm', vmin = vmin, vmax = vmax)
plt.show()

iband01.close()
gband01.close()
rband01.close()

Basic Analysis for M81

• Our first model does a very bad job matching our galaxy specifically for M81 as the current model we are using for sersic profile assumes a galaxy to be circular, wheras m81 is very clearly elliptical.
• When fit to an elliptical model, we obtain similar results to the other fits obtained using a one component-model, with our retrieved sersic index of 3.02 being reasonable for a bulge, but not capturing the broader spread of light within the larger disc structure.
• The two component model in this case gives slightly better results, with a sersic index of 1.75 for the disc, and strangely, .98 for the bulge.
• Fixing the index for the bulge yields even worse results, with an index of 3.11 for the disk, which is very unphysical for a spiral galaxy. Visually, the model seems to do a very poor job of capturing the outer halo of light from the galaxy. I hypothesize that this is because the inner disc is so prominent that even a two component model still does not capture the outer halo of light very well.

Model Fitting for M89

iband1 = fits.open("/Users/jackcolvin//m89_i-band_120.0s_bin1_250603_045958_jackcolvin_seo_0_FCAL.fits")
gband1 = fits.open("/Users/jackcolvin//m89_g-band_120.0s_bin1_250603_045721_jackcolvin_seo_0_FCAL.fits")
rband1 = fits.open("/Users/jackcolvin//m89_r-band_120.0s_bin1_250603_045449_jackcolvin_seo_0_FCAL.fits")

from astropy.visualization.lupton_rgb import make_lupton_rgb
from astropy.visualization import LogStretch


r1 = iband1[0].data
g1 = rband1[0].data
b1 = gband1[0].data

#lupton_rgb doesn't support scaling
#so that must be done manually
zscale = ZScaleInterval()
r_scaled1 = zscale(r1)
g_scaled1 = zscale(g1)
b_scaled1 = zscale(b1)



image = make_lupton_rgb(r_scaled1, g_scaled1, b_scaled1, Q = 1, stretch = .9, minimum = .3)
plt.imshow(image, origin = 'lower')
plt.show()

#Now we must find a best fit sersic profile
#and see if our best fit index fits our prediction 
#based on whether this galaxy is spiral or elliptical. 
#also check to see how different the function is across different bands
from scipy.optimize import least_squares
import matplotlib.pyplot as plt

scale_factor = 10000 #rescale calibrated fluxes as they are so small
cutout = Cutout2D(r1, (1050, 1050), (600, 600))
data0 = cutout.data
data = (data0-np.median(r1))*scale_factor
interval = ZScaleInterval()
vmin,vmax = interval.get_limits(data)
plt.imshow(data, origin = 'lower', cmap = 'gray', vmin = vmin, vmax = vmax )


y = np.zeros(shape = data.shape, dtype = int)
x = np.zeros(shape = data.shape, dtype = int)
for i in range(len(y[0])):
    y[i] += len(y[0]) - (i+1)
    x[:,i] += i

p0 = [np.max(data), 10, 4, 300, 300] #guess the expected sersic index for a spiral galaxy
result = least_squares(residuals, p0, args=(x, y, data), bounds=([0, 1, 0.1, 0, 0], [np.inf, np.inf, 10, 500, 500]))

I_e_fit, R_e_fit, n_fit, x0_fit, y0_fit = result.x
print(f"Fit: I_e = {I_e_fit:.5f}, R_e = {R_e_fit:.2f}, n = {n_fit:.2f}, x0 = {x0_fit:.2f}, y0 = {y0_fit:.2f}")

Fit: I_e = 1.13535, R_e = 17.60, n = 1.50, x0 = 266.66, y0 = 318.41

model_image = sersic_2d(x, y, *result.x)

import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm', vmin = vmin, vmax = vmax)
plt.show()

iband1.close()
gband1.close()
rband1.close()

Basic Analysis for M89

• There were much less considerations when fitting for elliptical galaxies than spiral, as for eliptical galaxies, the light is almost all concentrated within the central bulge, and there is not a prominent outer halo nor spiral arms to complicate the model. Therefore, I decided that a simple 1 component model would be sufficient to describe their structure.
• For M89, I got a sersic indec of 1.5, which is a significant understimation of the expected value of n~4 for elliptical galaxies.
• There are several physical reasons I noticed for this however. First of all, SEO has a fairly limited resolution, whcih causes light from point-like, dense sources of light to be spread out over a large area, which could account for the underestimation we see here
• Furthermore, the residual plot shows a similar trend, with the center of the model being an underestimation of the light intensity, and the outside being an overestimation. Therefore, it almost seems visually that a sharper model, and thus one with a larger sersic index, would be a better fit, despite what the residual fitting might return. This is certainly something to be investigated further as i continue to work on this project.

Model Fitting for m60

iband2 = fits.open("/Users/jackcolvin//m60_i-band_120.0s_bin1_250603_052218_jackcolvin_seo_0_FCAL.fits")
gband2 = fits.open("/Users/jackcolvin//m60_g-band_120.0s_bin1_250603_051942_jackcolvin_seo_0_FCAL.fits")
rband2 = fits.open("/Users/jackcolvin//m60_r-band_120.0s_bin1_250603_051649_jackcolvin_seo_0_FCAL.fits")

r2 = iband2[0].data
g2 = rband2[0].data
b2 = gband2[0].data

#lupton_rgb doesn't support scaling
#so that must be done manually
zscale = ZScaleInterval()
r_scaled2 = zscale(r2)
g_scaled2 = zscale(g2)
b_scaled2 = zscale(b2)



image = make_lupton_rgb(r_scaled2, .9*g_scaled2, b_scaled2, Q = 1, stretch = .9, minimum = .25)
plt.imshow(image, origin = 'lower')
plt.show()

#Now we must find a best fit sersic profile
#and see if our best fit index fits our prediction 
#based on whether this galaxy is spiral or elliptical. 
#also check to see how different the function is across different bands
from scipy.optimize import least_squares
import matplotlib.pyplot as plt

scale_factor = 10000 #rescale calibrated fluxes as they are so small
cutout = Cutout2D(r2, (1050, 1050), (600, 600))
data0 = cutout.data
data = (data0-np.median(r2))*scale_factor
interval = ZScaleInterval()
vmin,vmax = interval.get_limits(data)
plt.imshow(data, origin = 'lower', cmap = 'gray', vmin = vmin, vmax = vmax )


y = np.zeros(shape = data.shape, dtype = int)
x = np.zeros(shape = data.shape, dtype = int)
for i in range(len(y[0])):
    y[i] += len(y[0]) - (i+1)
    x[:,i] += i

p0 = [np.max(data), 10, 4, 300, 300] #guess the expected sersic index for a spiral galaxy
result = least_squares(residuals, p0, args=(x, y, data), bounds=([0, 1, 0.1, 0, 0], [np.inf, np.inf, 10, 500, 500]))

I_e_fit, R_e_fit, n_fit, x0_fit, y0_fit = result.x
print(f"Fit: I_e = {I_e_fit:.5f}, R_e = {R_e_fit:.2f}, n = {n_fit:.2f}, x0 = {x0_fit:.2f}, y0 = {y0_fit:.2f}")

Fit: I_e = 0.47759, R_e = 42.67, n = 2.00, x0 = 264.74, y0 = 325.45

model_image = sersic_2d(x, y, *result.x)

import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm', vmin = vmin, vmax = vmax)
plt.show()

#again, it seems like the fit is not capturing the true steepness of the real curve
#in a very similar way to the first elliptical. 

iband2.close()
gband2.close()
rband2.close()

Basic Analysis for M60

• For M60, I got a sersic indec of 2, which, while slightly better than my value for m89, still represents a significant deviation from the expectation for an elliptical galaxy.
• We again see the same pattern, where our residual plot suggests that a sharper model with a higher sersic profile might actually represent a better fit if some outluer points were removed.
• This galaxy is also slightly eliptical, so considering using an elliptical model instead of a circular one may be a good idea when I move forward with this project.

Model Fitting for M49

iband3 = fits.open("/Users/jackcolvin//m49_i-band_120.0s_bin1_250603_051201_jackcolvin_seo_0_FCAL.fits")
gband3 = fits.open("/Users/jackcolvin//m49_g-band_120.0s_bin1_250603_050926_jackcolvin_seo_0_FCAL.fits")
rband3 = fits.open("/Users/jackcolvin//m49_r-band_120.0s_bin1_250603_050648_jackcolvin_seo_0_FCAL.fits")

r3 = iband3[0].data
g3 = rband3[0].data
b3 = gband3[0].data

#lupton_rgb doesn't support scaling
#so that must be done manually
zscale = ZScaleInterval()
r_scaled3 = zscale(r3)
g_scaled3 = zscale(g3)
b_scaled3 = zscale(b3)



image = make_lupton_rgb(r_scaled3, .93*g_scaled3, b_scaled3, Q = .5, stretch = .9, minimum = .4)
plt.imshow(image, origin = 'lower')
plt.show()

#Now we must find a best fit sersic profile
#and see if our best fit index fits our prediction 
#based on whether this galaxy is spiral or elliptical. 
#also check to see how different the function is across different bands
from scipy.optimize import least_squares
import matplotlib.pyplot as plt

scale_factor = 10000 #rescale calibrated fluxes as they are so small
cutout = Cutout2D(r3, (1050, 1050), (600, 600))
data0 = cutout.data
data = (data0-np.median(r3))*scale_factor
interval = ZScaleInterval()
vmin,vmax = interval.get_limits(data)
plt.imshow(data, origin = 'lower', cmap = 'gray', vmin = vmin, vmax = vmax )


y = np.zeros(shape = data.shape, dtype = int)
x = np.zeros(shape = data.shape, dtype = int)
for i in range(len(y[0])):
    y[i] += len(y[0]) - (i+1)
    x[:,i] += i

p0 = [np.max(data), 10, 4, 300, 300] #guess the expected sersic index for a spiral galaxy
result = least_squares(residuals, p0, args=(x, y, data), bounds=([0, 1, 0.1, 0, 0], [np.inf, np.inf, 10, 500, 500]))

I_e_fit, R_e_fit, n_fit, x0_fit, y0_fit = result.x
print(f"Fit: I_e = {I_e_fit:.5f}, R_e = {R_e_fit:.2f}, n = {n_fit:.2f}, x0 = {x0_fit:.2f}, y0 = {y0_fit:.2f}")

Fit: I_e = 0.34682, R_e = 54.06, n = 2.09, x0 = 255.71, y0 = 309.56

model_image = sersic_2d(x, y, *result.x)

import matplotlib.pyplot as plt
plt.subplot(1, 3, 1)
plt.title('Data')
plt.imshow(data, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 2)
plt.title('Model')
plt.imshow(model_image, origin='lower', cmap='inferno', vmin = vmin, vmax = vmax)

plt.subplot(1, 3, 3)
plt.title('Residuals')
plt.imshow(data - model_image, origin='lower', cmap='coolwarm', vmin = vmin, vmax = vmax)
plt.show()

#again, it seems like the fit is not capturing the true steepness of the real curve
#in a very similar way to the first elliptical. 

iband3.close()
gband3.close()
rband3.close()

Basic Analysis for M89

• For M49, I got a sersic indec of 2.09, which agrees with my values for m60 and m89, and is slightly closer, but still significantly deviated from, the expected value of 4.
• Furthermore, I again see the same trend in my residuals that suggests a sharper model may be a better fit.

3. Discussions and conclusions

While I was unable to succesfully reach the first goal of this project – matching the expected sersic frofiles of n~1 for spiral galaxies and n~4 for ellipticals, I did still have some successes. After using a two component model for disk galaxies, I was able to measure a value of ~1 for 2/3 of my spiral galaxies, and n>2 for all three. Furthermore, I did measure elliptical galaxies to have a more dense light distribution than spirals, with all three eliptical galaxies having n~2.

I had to deal with significant sources of error throughput my analysis for both eliptical and spiral galaxies. In both cases, the tendency of SEO to spread light across a larger number of pixels due to its limited resolution (3-4 arcseconds) certainly had an impact on my fits. For spiral galaxies, the different light profiles of both the bulge and the spiral arms was a large source of errors that biased by values, and I was only able to partially mitigate those error source by a combination of masking, smoothing, and making a more complex, two-component model. For elliptical galaxies, The limited angular resolution of SEO likely had an even larger effect, as their light is much more similar to a point source, and thus requires high angular resolution to accurately measure. Furthermore, I observed a systematic underapproximation of sersic indeces in all three galaxies, evidenced by residual plots that sugested that a sharper model may actually be a better fit, at least for the central core of the galaxy.

As I continue to work on this project, I have a couple key considerations to keep in mind. First of all, I want to continue tweaking my two-component spiral galaxy model to try to retrieve a model that can simulataneously give me accurate values for both my bulge and my disk. Furthermore, I wasnt to work more on smoothing to lessen the effect that the spiral arms may be having on my sersic fits. For the eliptical galaxies, I want to try to contrain my other parameters apart from sersic index more tightly to see if i can retrieve a model that qualitatively reduced my residuals more effectively. Furthermore, for both cases, I want to both take more SEO eposures, and maybe also try some higher reolution images from the Legacy or SDSS surveys, so try to isolate whether my errors are a result of SEO's poor resolution, or some aspect of my specific analytical techniques.